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The electrostatic potential in a superconductor is studied. To this end Bardeen's extension of the 
Ginzburg-Landau theory to low temperatures is used to derive three Ginzburg-Landau equations 
- the Maxwell equation for the vector potential, the Schrodinger equation for the wave function 
and the Poisson equation for the electrostatic potential. The electrostatic and the thermodynamic 
potential compensate each other to a great extent resulting into an effective potential acting on 
the superconducting condensate. For the Abrikosov vortex lattice in Niobium, numerical solutions 
are presented and the different contributions to the electrostatic potential and the related charge 
distribution are discussed. 



I. INTRODUCTION 

Even in equilibrium, any inhomogeneous conductor has 
internal electric fields which keep its charge distribution 
close to local neutrality. The superconductor is not an 
exception. While the electrochemical potential is con- 
stant, the local chemical potential varies in general with 
any gradient in the system. A distinct property of the 
superconductor is that in equilibrium there can be an 
inhomogeneity due to the diamagnetic electric current. 

The electric field in a superconductor with a statio: 
ary current has been discussed already in 1937 by BoppI 
From the hydrodynamic description of a charged liquid, 
Bopp has concluded that the inertial and Lorentz force 
created by the current are balanced by the Coulomb 
force. The corresponding electrostatic potential has the 
form of a Bernoulli potentialH. 

If the Lorentz force dominates, the Bernoulli poten- 
tial can also be considered as Hall effect. While it was 
clear that there has to be a Hall voltage which passes the 
Lorentz force from electrons to the lattice, its measure- 
ments by contacts in standard Hall setups did not show 
any. It was understoodu that by contacts one observes 
differences in the electrochemical (not electrostatic) po- 
tential but this potential is constant in equilibrium. 

With the aim to distinguish the electrostatic potential 
from -the electrochemical one, as late as 1968, Bok and 
KleinQ have used the Kelvin capacitive coupling proposed 
by Huntl3 and have observed first the Bernoulli potential 
on the surface of a superconductor. Similar measure- 
ments have been performed by Brown aiid. MorriaJij or 
more recently by Chiang and ShevchenkdZllj. 

Even a perfect surface establishes itself a very stmng 
defect which essentially modifies the electric field.El It 
is desirable to observe the internal electric field directly 
in the bulk of a superconductor. A new experiment in 
this direction has been performed recently by Kumagai 



et aZ0 who have measured the electric field in a type- 
II superconductor in mixed state by nuclear quadrupole 
resonance. 

Another consequence of the electric fieldj-Ln the bulk 
is a charge of the vortex core. Blatter et aZEJ have pro- 
posed an experiment by which the vortex charge can be 
accessed. Such measurement, however, is still to be per- 
formed. It is also speculated that the vortex charge af- 
fects the motion of vortices and, thus plays a role in the 
sign reversal of the Hall regimell3. Since the theory of the 
anomalous Hall voltage is still open, one cannot conclude 
about the core charge from this effect. 

In this paper we derive a phenomenological theory of 
the Ginzburg-Landau (GL) type which allows one to eval- 
uate the electric field in the bulk of superconductors at 
low temperatures. A brief presentation of this theory has 
been already published in Ref. |l^. Here we present de- 
tails and show how to handle numerically this theory for 
the Abrikosov lattice of vortices. The electrostatic poten- 
tial in the vortex lattice is shown for a selected temper- 
ature and the contribution of the electric field to forces 
acting on the condensate is discussed. Throughout the 
paper we use the language of the two-fluid model. The 
fluid of superconducting electrons is called condensate 
while electrons mean normal electrons. 

In the next section we re view t heoretical approaches 
to the electric field. In Sec. [II A we introduce the free 
energy which includes the condensation energy of Gorter 
and Casimir, the kinetic energy of Ginzburg and Lan- 
dau, and the standard electromagnetic energy. Sec. [II B 
presents the essential part of our approach. We use the 
variational principle to derive three GL equations: the 
Maxwell equation for the magnetic field, the Schrodinger 
equation for the wave function, and the Poisson equation 
for the electrostatic potential in the bulk of supercon- 
ductors. In Sec. the hydrodynamic picture is used 
to link the presented theory with the former approaches 



reviewed in Sec. ||. In Sec. ^ we discuss magnetic prop- 
erties of the Abrikosov vortex lattice as a function of the 
temperature. In Sec. ^ we compare the electrostatic 
potential with other potentials acting on the condensate. 
We also present the charge distribution and show that its 
amplitude is very small what allows one to employ a con- 
venient quasi-neutral approximation. Sec. VI] presents 
the conclusions. In Appendix A we estimate the mate- 
rial parameters for Niobium using the McMillan formula 
and empirical rules established from chemical trends. 



II. HISTORICAL REVIEW 

The electric field in superconductors has been studied 
since the discovery of superconductivity. Accordingly, 
various approaches to this problem can be found in the 
literature. We will briefly remind the progress in this 
field made mainly in late 1960's and early 1970's. 



VjVjAi = — [v X V X A]i + VjViAj 



(4) 



In the first term of one can recognize the Lorentz 
force, ev X V X A = ev X B. In the second term of we 
substitute A by the velocity from the London condition, 
eVjViAj = — mvjViVj — —Vi^mv^. 

The time derivative of the London condition then reads 

' (5) 



mv = e(E + V X B) + V ( e</j + -mv 



This equation can be compared with the Newton equa- 
tion (0) giving the electrostatic potential as 



1 9 

Vetp = Fs - V-mw . 



2. Bernoulli potential 



(6) 



A. Bernoulli potential 

The Bernoulli potential for superconductors has been 
first derived by BoppEI. Here we follow the later approach 
of Londona. The condensate has to obey two equations 
of motion. First, it is the London condition. 



= -eA, 



(1) 



where v is the local velocity of the condensate and A is 
the vector potential. Second, it is the Newton equation 



mv = e(E + V X B) 



(2) 



where the first term is the Lorentz force with the electric 
field, E = —dA/dt — V(/?, and the magnetic field, B = 
V X A. The additional force has been treated by 
different authors within rather different approximations. 

Since the motion of the condensate is fully determined 
by the London condition, one can use the Newton equa- 
tion to determine the force acting on the condensate. 
Once the additional force will be specified, this proce- 
dure allows one to identify the electrostatic potential ip. 



1. Time derivative of the London condition 

To bring the London condition into a form which can 
be easily compared with the Newton equation, we take 
the total time derivative, d/dt = d/dt -f- (vV), of the 
London condition (|^), 



dA 

mv — —e— e(vV)A. 

dt ^ ' 



(3) 



The first term we express via the electric field, —dA/dt = 
E -|- V(p. For the second term we use a vector identity 
which in components reads 



London assumed that the motion of the condensate is 
controlled by the Lorentz force only. In this approxima- 
tion, there is no additional force, 



0. 



(7) 



From (||) thus follows the electrostatic potential of the 
Bernoulli type. 
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-—mv . 
2 



(8) 



3. Quasiparticle screening 

In 1964 van Vijfeijken and Staas0 have extended the 
Bernoulli potential to finite temperatures using the two 
fluid model. When flowing, normal electrons dissipate 
energy. Therefore, in the stationary case they have to 
stay at rest in spite of the presence of an electric field. 
These authors have introduced an unspecified force, 



eV(/3, 



(9) 



acting on electrons to keep them at rest, F„ -I- eE = 0. 
This force is assumed to result from the interaction be- 
tween the electrons and the condensate. Accordingly, 
there has to be a reaction force F^ acting on the con- 
densate so that the Newton law of action and reaction is 
fulfilled. 



UsFs = 0, 



(10) 



where n„ and Hs are densities of electrons and conden- 
sate. From (Q) and (10) one finds the additional force, 



Fs = eVip, 



(11) 



2 



and from m) follows the electrostatic potential 



ris 1 
~n2' 



(12) 



This is the Bernoulli potential (j^) reduced by the share 
of the condensate on the total density, n — Un + Us- 

The reduction of the Bernoulli potential has become 
known as "screening by normal electrons" or "quasipar- 
ticle screening" . The quasiparticle screening, however, 
has to be distinguished from the Thomas-Fermi screen- 
ing present in all metals including superconductors. 



4-. Thomas-Fermi screening 

In superconductors, the screening is the same as in nor- 
mal metals. Starting from the time-dependent Ginzburg- 
Landau theory, Jakeman and Pikell3 have derived the 
Poisson equation for the electric field with the reduced 
Bernoulli potential as the driving term. 



ns_l 
~2' 



(13) 



Currents change typically on the scale of the London 
penetration depth or the GL coherence length, which 
are much larger than the Thomas-Fermi screening length 
Xtf- The electrostatic potential (p thus can be treated 
in the limit of strong screening, Xtf ~* 0, and from ( p^ 
one recovers (f2h. 



5. Thermodynamic potential 

Already in 1949 Sorokin0 has followed the hydrody- 
namic approach of Bopp assuming an unspecified free 
energy, 



:Fs = / drf, 



(14) 



responsible for the superconducting transition. Here fs 
is the density of free energy and dr denotes integration 
over the sample volume. The free energy leads to a ther- 
modynamic potential. 



(15) 



Shs dris ' 
which yields the additional force 

F, = -Vws. (16) 
According to (^) the Bernoulli potential is modified as 
1 o 



ef = ——mv — Wg 



(17) 



The quasiparticle screening is one of the contributions 
that result from the thermodynamic potential. There are 



also other contributions which can provide information 
about the pairing mechanism. 

Unfortunately, London has disregarded the thermody- 
namic potential in his bookH as unknown and unimpor- 
tant. His objection was correct at that time since the 
first reliable thermodynamic potential has been derived 
eight years later by Bardeen, Cooper and SchrieffertZl. 
On the othe&jliapd, the two- fluid free energy of Gorter 
and Casimiilia~E3 known from 1934, could be used within 
Sorokin's approach to provide at least qualitative results. 
Our approach follows Sorokin, except that we use an ex- 
plicit thermodynamic potential of Gorter and Casimir 
and a non-local kinetic energy. 



6. Non-local corrections 

As shown in Ref. |l|, London's approach can be mod- 
ified towards strongly inhomogeneous systems using the 
Schrodinger equation for a Cooper pair. 



2m 



0, 



(18) 



instead of the Newton equation (g). Here we have also 
included the thermodynamic potential Ws neglected in 
Ref. m. 

From ( [T^ ) follows directly a quantum modification of 
the Bernoulli potential. 



^ -0 2™* 



2w... 



(19) 



In the quasi-classical approximation, (—ihV — e*A)?/' — 
m*vip, this formula reduces to potential ( |l7| ) derived by 
Sorokin. 

To obtain the actual value of the potential, the wave 
function -0 is identified with the GL wave function and 
solved from the GL equation. Accordingly, the Cooperon 
mass and charge, m* = 2m and e* — 2e, appear in the 
Schrodinger equation (O). 



B. Thermodynamic correction 

Rickayzenil proposed a thermodynamic approach to 
the electric field. He assumes a quadratic dependence of 
the free energy on the velocity, what limits his study to 
weak currents. For systems with a parabolic band, the 
increase of the free energy due to the current equals the 
kinetic energy of the condensate. 



/kin 
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Us-mv 
2 



(20) 



The electrochemical potential, v = Ep -\- i^kin + e^j 
is constant in the whole system, therefore et/? = — J^kin- 
Since v = df jdn^ the velocity variation of the local chem- 
ical potential is i^Mn = dfua/dn. Accordingly, ±he elec- 
trostatic potential induced by the current readstJ 



3 



dus 1 



(21) 



Expression (^) generalizes (^2|). From the phe- 
nomenological density of the condensate, 



1 



follows 



eip — 



—mv 
n 2 



n„ 9 In Tc 1 , 
-—mv . 



n 9 Inn 2 



(22) 



(23) 



The first term is the reduced Bernoulli potential (p^), 
the second is a thermodynamic correction. According to 
(^2|), the first term of ( p3|) depends on the temperature 
as 1 - T'^/T^ while the second one goes as T^/T^. At 
higher temperatures the second term dominates. 



1. BCS estimate 

The density dependence of reflects the pairing mech- 
anism, pits magnitude can be estimated from the BCS 
relationliZi, 



(24) 



where T) is the single-spin density of states, is the cut- 
off frequency usually approximated by the Debye tem- 
perature, hiOD ~ kBOn, and V is the BCS interaction. 
Assuming that 6d and V do not depend on the density, 
one finds 



dlnTc 
9 In 71 



dV 



1 



dlnnV^V 



' a Inn ^n' 



(25) 



It remains to estimate the derivative of the density of 
states. For systems with a parabolic band the density 
of states is proportional to the Fermi momentum. In I? cx 
kp, while the density of electrons isn oc kp. Accordingly, 
dlnV/dhin « 1/3. For Niobium we have a very similar 
value In 13/(9 Inn = 0.32, see Tab. 1 in Appendix A. 
With the BCS estimate (|5|), the electrostatic potential 
23|) reads 



z V n n 6 In 



(26) 



For conventional superconductors, O^/Tc is of the order 
of few tens, therefore the thermodynamic correction is 
the dominant contribution for approximately T > ^T^. 

For Niobium the BCS formula ( p6| ) overestimates the 
thermodynamic correction. The approximate factor from 
(|§ is (4/3)ln(6lD/'rc) = 4.5 while the full factor from 
(|3| gives -4(91nTc/91nn) = 3.0, see (|§) and Ap- 
pendix A. 



2. BCS microscopic theory 

Within the BCS theory, the-filectric field-has been stud.- 
ied by Adkins and Waldrarro, RickayzenE3 and Hongtil. 
In all these studies, materials with a general band struc- 
ture have been addressed. For the sake of simplicity we 
discuss onljfr-the parabolic band, for which the BCS the- 
ory yieldsEScJ 



(27) 



Here Aq = UdhsTc is the gap at T = and A is the 
actual local value of the gap. 

Since the electric current locally depresses the gap. 



A = A 



oq 



A' with A' 



the potential {UTh in- 



cludes the contribution of Bernoulli type, up. = (pcq + ip' 
with (p' (X —v^. As shown by RickayzenEj, can be 
rearranged into the thermodynamic correction of (E9). 



C. Aims of the present approach 

In this paper we discuss the Ginzburg-Landau theory 
modified in two directions. First, following Bardeen we 
use its extension to low temperatures. Second, we in- 
clude the electrostatic potential. We focus on the bulk of 
superconductors, i.e., on regions which are far from the 
surface on the scale of Thomas-Fermi screening length. 

Starting from the free energy, we derive the Poisson 
equation along with the Maxwell equation for the vector 
potential and the equation of the Schrodinger type for 
the wave function. The presented theory yields 

• non-local Bernoulli potential, 

• quasiparticle screening, 

• thermodynamic corrections, 

• thermoelectric field of normal metal at T = Tc, 

• Thomas-Fermi screening. 

Our approach parallels the original study of Sorokin, 
however, we use the explicit phenomenological free en- 
ergy proposed by Bardeen. It combines the Ginzburg- 
Landau (GL) theory with the Gorter-Casimir free energy. 
Naturally, this theory is only approximate. Its major ad- 
vantage is its transparency and a simple implementation 
scheme. 



III. EXTENDED GINZBURG-LANDAU THEORY 



A. Free energy 



BardeenBE^ has extended the GL theonf|ffl by the 
use of the Gorter-Casimir two-fiuid mode£a~E3 so as to 



4 



apply to all temperatures. We briefly recall the Gorter- 
Casimir model and introduce other components of the 
free energy. 



1. Condensation energy of two-fluid model 

Gorter and Casimir assumed that the superconducting 
state is characterized by an order parameter w which is 
zero in the normal state and unity at zero temperature. 
They have modified the normal state density of free en- 
ergy as 



fs^U-e. 



1 



con^ - -'-fT'^VT 



ru. 



(28) 



For w = 0, the free energy (|2^) equals the normal state 
free energy consisting of the internal energy U and the 
entropy term —^'yT'^. Sommerfeld's 7 is the linear coef- 
ficient of the specific heat. In the superconducting state, 
vu 0, two mechanisms are expected. First, the order- 
ing releases the condensation energy £conW. Second, the 
ordering reduces the entropy by the factor VI ~ 

In equilibrium the free energy reaches its minimum. 
From STg/Sm = dfs/dm = follows that the equilib- 
rium value of the order parameter is a solution of 



4^/^ 



(29) 



At the critical temperature the ordering vanishes, m — 0, 
therefore 



1 2 

£con = ^T^c 



From (g9|) with (|30|) follows 

VJ = 1 7, 



(30) 



(31) 



which agrees with the observed temperature dependence 
of the condensate density (p^). Accordingly, one can 
identify the order parameter w with the fraction of the 
condensate on the total density of electrons. 



n 



(32) 



2. Kinetic energy 

An electric current contributes to the free energy by 
the kinetic energy of the condensat£L.JIhe kinetic energy 
proposed by Ginzburg and LandauEZIca reads 



dr- 



2m* 



-iW- e*A) 



(33) 



(34) 



We use the isotropic effective mass for simplicity, the 
anisotropic case will be discussed in the next paper. 

In the Gorter-Casimir free energy we thus substitute 
the order parameter by 



(35) 



3. Electromagnetic energy 



The density of free energy has four components. 



(36) 



The free energy JF, given by the volume integral (|T^) of 
the free energy density (|2^), we shall call the condensa- 
tion energy according to its most important part. The 
kinetic energy, JFkin, is given by the GL expression (|33|). 
The Coulomb interaction reads 



1 



1 



J^c = l [[ drdr' ^ , 

2 JJ Aire |r - r 



-p{r)p[r'), 



(37) 



where p = e*\'ip\^ + eun + putt is the charge density. The 
Coulomb interaction also determines the electrostatic po- 
tential, by 



^(r) 



dv' 



1 



1 



47re r — r' 



(38) 



or in its differential form by the Poisson equation 

- eVV = P- (39) 
Finally, the Helmholtz magnetic free energy readsS, 



A[ 



dr-^(B-Ba)2, 



(40) 



where is the applied magnetic field. 

The total free energy is a functional of the wave func- 
tion, the vector potential and the (normal) electron den- 
sity, J^[V', A, n„]. The other physical quantities like B, 
n, Us, p ot: VJ are subsidiary and have to be understood 
as functions of the independent variables t/i, A and n„. 

We note that much more sophisticated approximations 
of the free energy J^s have been developed from the 

BCS theory and Eliashberg's theory already in 1960's, see 
e.g. Ref. In principle, one can start from any of these 
approximations. Since our prime interest is in the elec- 
trostatic potential and the related charge distribution, 
we prefer to use the simple approximation of Bardeen. 



with the wave function normalized as 



5 



B. Ginzburg-Landau equations of motions 

In equilibrium, the system stays in the state with min- 
imum free energy. Accordingly, the variations of T with 
respect to the vector potential A, the wave function '\\) 
and the electron density n„ have to vanish. 

During the variation procedure, the two-point func- 
tion Tc and the one-point functions (all the others) are 
treated differently. The local contributions are given by 
the corresponding densities, /„[/(r), V/(r)], witb—Fo, ~ 
J drfa, and their variation is of Lagrange's formO, 



SI dl dVl ' 



(41) 



Here a represents the subscripts s, kin and M while I 
stands for A, or n„. 

The variation of the Coulomb energy with respect to 
"0 or n„ can be expressed by the variation with respect 
to the density of charge p which reads 



6p{r) 



= / dr' 



1 



1 



Aire r — r' 



(42) 



According to ( p^ ) we can abbreviate this variation as 



Sp 



(p. 



(43) 



The variation with respect to ip leads to the equation 
of the Schrodinger type,[| 



1 



2m 

The effective potential, 
S 



X 



{Tc + Ts) = e*ip- 



dfs 



(46) 



(47) 



covers all forces acting on Cooper pairs. 

3. Diffusion of normal electrons 

From the variation with respect to the electron density, 
8T I bun = 0, one finds that the sum of all potentials 
acting on the normal electrons has to vanish, i.e., 



5Ts 
Snn 



dfs 

dUn 



(48) 



This condition parallels Eq. (^ of van Vijfeijken and 
Staas. 

The set of equations ( p^^^ ) is closed. Its particular 
form is given by the condensation energy fg. Below we 
evaluate derivatives of the condensation energy within 
the Gorter-Casimir approximation (ESh. 



1. Maxwell equation 



4- Effective potential acting on Cooper pairs 



The vector potential A appears in the kinetic energy 
^kiii) and its gradients enter the magnetic free energy via 
B = V X A. From the condition of minimum with respect 
to A, one recovers the Maxwell equation. 



V X V X A = ^oj, 



(44) 



where tiie current j is given by the quantum-mechanical 
formulaES, 



j = —Reij{-ihV-e*A)i^, 
m 



(45) 



known as the second GL equation. Here ip denotes the 
complex conjugate of ip. 



2. Schrodinger equation 



To describe the motion of the condensate given by the 
effective potential x, we have to evaluate the electrostatic 
potential. This will be done in the spirit of van Vijfeijken 
and Staas using the equation for normal electrons (HqV 

Since e* — 2e, the effective potential results from ([48|) 
and (1^ as 



X 



dfs 



,5^ 

dUn 



(49) 



The combination of the partial derivatives excludes a 
contribution of functions which depend exclusively on the 
total density, dn/d\ip\'^ — 29n/9n„ = 0. Accordingly, 
derivatives of density-dependent material parameters {U, 
^con, 7 and Tc) do not contribute to the potential %. From 
ll) and <M) thus follows 



The wave function ^ enters the free energy via the or- 
der parameter, the Coulomb interaction via the charge 
density, and the kinetic energy J^kin, where also the gra- 
dients of ijj appear. ThepSiariation parallels the derivation 
of the first GL equatioECJ, for details see Ref. E9^. 



* The free energy is a real function which depends on the 
complex function ^ and its conjugate tp. We express absolute 
values as products, |?/>p — tptp, and take 5'tp and Sip as inde- 
pendent perturbations, 5T = {5T /5il))Sil) + {ST / 5'ip)5'ip. Since 
5J-/S^ = 5T/5%1), both variations vanish at the same time. 
The variational condition 5T /S^p — yields the Schrodinger 
equation. 
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X = -2 + -TT- 

n 2n 



(50) 



With the potential (pO|), the Schrodinger equation ( ^6| ) 
is identi£:a,l_to the extended GL equation proposed by 
BardeenBEl. 

Close to Tc, the potential x approaches the quadratic 
form of Ginzburg and Landau, 



with 



X 



n 



= — 



(51) 



We note that the effective potential x depends on the 
density n via 7 and Tc- In principle, one has to iterate the 
GL equation together with relations for the density n. In 
practice, deviations of the density from its crystal value 
are very small, |n + putt/e-l ^ n, and the approximation 
en « — piatt is well justified when one solves for ip and A. 



5. Poisson equation with screening 



Now we rearrange (^8|) into the form of the Poisson 
equation with the Thomas-Fermi screening. The varia- 
tion of the free energy in (Eq) reads 



dUn 



dU_ 
dn 



d 

drir, 



The derivative of the internal energy. 



dn 



Ep, 



(53) 



(54) 



is the Fermi energy at the normal ground state of total 
electron density n. In the presence of the electrostatic 
potential, the density n differs from its crystal value, 
no = — piatt/e, by the density perturbation p/e. The 
Fermi energy thus depends on the charge density as 



Ep — Ep 



dEp p 
dn e ' 



(55) 



where Ep is the crystal value. As usual in the theory of 
superconductivity, we associate the crystal Fermi energy 
with the origin of the energy scale, Ep — 0. 

The density dependence of the local Fermi energy de- 
termines the screening. The density derivative of the 
Fermi energy is the inverse density of states. 



dEp 
dn 



1 

2D' 



(56) 



Using the Poisson equation (|39|) to evaluate the charge 
from the electrostatic potential, p = —e^'^ip, one can 
express the Fermi energy as 



E, 



_y2 V72 



'\pp 



(57) 



with the Thomas-Fermi screening length \^p = 
e/{2Ve^). 

By a substitution of the Fermi ener gy ( |57| ) in the sta- 
bility condition for normal electrons (Eq), we arrive at 
the screened Poisson equation. 



dnn 



(58) 



The right hand side of ( p8| ) is readily evaluated from 
the assumption that the condensation energy econ and 
the Sommerfeld 7 depend only on the total density 



and from the explicit form of the or- 



der parameter, w = 21-01^/(21-012 _|_ ^^-^^ giving 



Arpp 



X- 



|0|2 9£con2|V|- 



dn 



2101^ 



(59) 



In the language of Jakeman and Pike, the Poisson equa- 
tion ( |59| ) is called third GL equation. 

The first term on the rhs of (^^ is the non-local 
Bernoulli potential with quasiparticle screening. This 
can be seen if we multiply the GL equation ( |46| ) by ip 
which yields 



X- 



2m*n 



0(-i?iV-e*A)^V- 



(60) 



In the classical approximation of the kinetic energy, 
{l/2m*)ilj{-ihV - e*Ay'tp « \m* v'^\'ip\'^ , one finds that 
the first term of (|5^) is the screened Bernoulli potential 
of van Vijfeijken and Staas, xl^P/n ~ —{ns/n)^nnv'^. 

The second and third terms of ( p9| ) are non-linear gen- 
eralizations of the thermodynamic correction by Rick- 
ayzen. Note that the third term remains finite at the 
critical point, T Tc and. 101 — > 0, yielding the normal 
state thermoelectric fieldca. 

The set of GL equations is closed. It consists of 
the Maxwell equation (|44| ) with the current ([45|) , the 
Schrodinger equation ([46|) with the potential (^0^ , and 
the screened Poisson equation ( |59| ) . Deviations from the 
local charge neutrality are given by the bare Poisson 
equation (|39|). 



IV. HYDRODYNAMIC PICTURE 



Within the thermodynamic approach of Sec. IIIB, the 
electrostatic potential ip is a function of the wave func- 
tion ip. This contrasts with the original derivations of the 
Bernoulli potential expressed in terms of the condensate 
velocity v. To make the link with the original approaches 
mentioned in Sec. ||, in this section we reformulate the 
above thermodynamic theory in the hydrodynamic pic- 
ture. 
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The hydrodynamic picture is readily obtained writing 
the wave function in terms of the condensate density (Q) 
and the phase 9, 



(61) 



A velocity defined via the current, j — en-jV, then reads 



V = — (hSie - e*A) 

m* 



A. On Sorokin's relation 



(62) 



In the representation (|6^) the Schrodinger equation 
1) reads[] 



— mv — — 

2 8to ^/nl 



frTs + + Wa = 0, 



(63) 



where we have used m* = 2m and x — 2e(^ + 2ws as it 
follows from (^^. The thermodynamic potential Wg is 
given by 



1 9/, 



2 9|V'p 



(64) 



which is equivalent to (|T^). With the non-local correc- 
tion neglected, \J'^^/nl k. 0, equation (|6^ ) turns into the 
Sorokin result (17). 

Naturally, the explicit evaluation of the electrostatic 
potential within the Sorokin approach parallels the non- 
local approach presented in the previous section. Since 
this approach is quite transparent, we show this proce- 
dure in detail. 



Relation (50) reads 



-con 

n 



An 



V n 



(65) 



Substituting ( |65[ ) into the local approximation of (|63[) one 
finds 



1 



-mv 



1 



4n Jl-^' 

V n 



(66) 



which yields the condensate density as a function of the 
local velocity. Provided that the profile of velocities in 
the system is known, from (^ and (|5^) one can directly 
evaluate the electrostatic potential. 



^Equation (^s]) is the energy- conserving integral of motion 
of the Newton-like form of the Schrodinger equation. This 
Newton-like equation itsptf may be found e.g. in The Feyn- 
man Lectures on PhysicsE3. 



B. On Rickayzen's result 

Now we recover Rickayzen's result (p3|). To this end we 
have to accept identical approximations. First we neglect 
the Thomas-Fermi screening, so that (M) reads 



eip = — {eip + Ws)^ 

n n on 



1 



n 



(67) 



Second, Rickayzen assumes a local relation between the 
velocity and the electrostatic potential. Accordingly, we 
neglect the gradient correction in (|6^), i.e., we use the 
Sorokin approximation, eip + Ws = — (1/2)tow^, with the 
help of which wc eliminate the Sorokin potential Ws from 



Third, following Rickayzen we take the limit of weak 
currents, nmv^/2 ^ £con -_7T^/4. Up to hnear orders 
in the kinetic energy, from (66) follows — n'^ + n'^ with 



rp4 



(68) 



so that from (^ 

<Peq + as 



results the electrostatic potential Lp 



n„ 1 



eip — 



—mv 
n 2 



dr 



dj 1 n'^T^ 



(69) 



Using (^) and ([2^) , one can rearrange expression ( |69| ) 



into the potential (E3[) derived by Rickayzen. 

In summary, the result of Sorokin corresponds to the 
local approximation of the presented approach. The po- 
tential of Bernoulli type derived by Rickayzen includes 
further approximations, in particular the limit of the 
weak electric current. We note that for systems with 
vortices the non-local approach is necessary since the 
"classical" kinetic energy mv"^ /2 diverges at the vortex 
center. This divergence is compensated by the non-local 
correction so that the "quantum" kinetic energy remains 
regular. 



V. MAGNETIC PROPERTIES OF THE 
ABRIKOSOV VORTEX LATTICE 

In this section we evaluate the wave function and the 
magnetic field for the Abrikosov vortex lattice in Nio- 
bium. Pure Niobium is close to the border between type- 
I and type-II superconductors since its GL parameter 
K = 0.78 is only shghtly above l/v^. However, the GL 
parameter can be increased up to about three by impu- 
rities. For simplicity we neglect the effect of impurities 
on material parameters other than the GL parameter n. 

As will be proven in the next section, deviations of the 
total density of electrons from its unperturbed value are 
very small, |p| <C piatt- We will neglect these deviations 
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and treat the material parameters in the approximation 
of quasi-neutraUty, j{n) « 7(no) etc. In this approxi- 
mation, the first and the second GL equations are inde- 
pendent of the third GL equation. Therefore, we shall 
ignore the electrostatic potential and related charge de- 
viation within this section. 



A. Dimensionless notation 

Our approach parallels Ref. [sj. In our calculations we 
shall use dimensionless quantities, 



t 



T 



A 

r 



Alo 



(70) 



Close to the critical temperature, t I, (these dimen- 
sionless variables reduce to the usual form.c3 

The thermodynamical critical field Be, the London 
penetration depth Aloii, and GL parameter k depend on 
the temperature as 

B,{t) = Bail-t''), 
Ao 



Alo„(<) = 



K(t) = Kq 



l+t2' 



Bc2{t) = V2kBc = 2Bq 



Kq 



(71) 



The asymptotic values of these quantities in terms of the 
parameters of the Gorter-Casimir model read 



Bo^T, 
Ao 



W7 
2 ' 



neh y jiQ 



(72) 



Finally, we introduce a dimensionless amplitude of the 
wave function and the dimensionless velocity. 



2m 



n(l-t4)^ 
Q = a - -V6I. 



(73) 



Our dimensionless notation is identical to Ref. |3J. 

The Schrodinger equation (^6|) with the effective po- 
tential (po|) in the dimensionless notation reads 



— v2tj + VV^ 
( 



1 -f2 



v/1 - (1 - t^)^ 



- 1 



(74) 



The terms on the left hand side result from the kinetic 
energy, the terms on the right hand side represent the 
potential. 

The Maxwell equation (|4^) with the current (|4^) reads 



V^Qb = -ujQ,A - t^Qb 



(75) 



The full quantum velocity is a sum of two terms, Q = 
+ Qf), where Qa is any model velocity field which 
covers the singular contributions at vortices. 



V X = b - $0 51 '^(f ~ ^) 



(76) 



R 



In our choice of Q,a appears the mean value of the mag- 
netic field in the superconductor. 



b = (b> = ^ / dfb. 



(77) 



therefore, for an ideally periodic vpstex lattice Q,a is 
given by the Abrikosov Bc2 solution.Eil The sum over the 
2D (5-functions represents the contributions of the nodes 
of the wave function in the vortex centers R to the quan- 
tum velocity. In this choice one has (V x Q^) = so that 
V X Qft = b — b describes the spatial modulation of the 
magnetic field due to the diamagnetic currents. 



B. Fourier representation 

For the periodic lattice of vortices, it is advantageous 
to express all functions by Fourier series. 



uj{r) = ^ aK(l — cosKf), 
5(f) ~b + ^ 5k cos Kf . 



(78) 
(79) 



We choose the direction of vortices along the ^-axis. The 
function 5(f) is the z component of the magnetic field b. 
Since the system is translationally invariant along z, the 
vectors f = (a;, y) and K (reciprocal lattice vectors) are 
two dimensional. 

The special choice, w oc (1 — cosKf), enforces nodes of 
the wave function, tt'(R) = 0, at the positions of the vor- 
tex centers, R = iixi+jx2, jy2) with i, j — 0, ±1, ±2, . . .. 
For the triangular lattice, two of the nearest neighbor 
vortices are at R = {xi,0) and R = (2^2, 2/2), where 
Xi = 2x2 and 2/2 = ^f?>X2- The distance between vor- 
tices, xi, is determined by the condition that each vortex 
contributes to the mean magnetic field by one elementary 
quantum of flux, Sh = Xiy2b = $0 = 2tt/k. The sums 
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in the Fourier representation run over nonzero discrete 
momenta, K = {2-k / S){iy2: jxi + 1x2). In this choice of 
the Fourier expansion the mean value of the amphtude 
of the wave function reads 



a; = (cj) = ^ CK- 
Since Q;, = Q — is periodic one may write 



Q(f) = QA(f) + &K^^sinKf 



(80) 



(81) 



with z X K = {Ky, —Kx) and z is the unit vector along 
the axis z. 



C. Simple iteration scheme 

Now we are ready to specify the iteration scheme for 
the Fourier components of the wave function and the 
quantum velocity. 

The Fourier representation of Eq. ( fz^ ) reads 



Ok 



((s -2uj + ujQ'^ + g) cosKf) 



A'2 



with 



1 



1 



The Fourier representation of Eq. ^73) reads 



5K 



2 {{ujb + il>{b - b) + p) cos Kf ) 



with 



p = (Vw X Q)z = Qx- Qyiz-- 

oy ox 



(82) 

(83) 
(84) 

(85) 
(86) 



Within a simple iteration scheme for given values of 
t, b, and k, one starts from the Abrikosov Bc2 solution 
or some other values of ok and 5k. In the step (a) one 
evaluates hk from ( ^2|) and upgrades uj{x,y) according 
to ([78|). In the step (b) one evaluates 5k from (p^) and 
upgrades b{x,y) and Q{x,y) from ( [79|) and (^1]). The 
iteration scheme (a), (b), (a), (b), . . . then leads to the 
periodic solution of Eqs. (Q) and (|75|). 



each use of equation (82). Here we show how to make 
this optimization within the Bardeen set of equations. 

Assume a change of the wave function Lu{x,y) which 
maintains its shape but modifies its amplitude, 



<^{x,y) := {l + c)uj{x,y), 



(87) 



i.e., the old value u; (the right hand side) obtained from 
( p2[ ) is replaced by a new value. 

The constant c has to be found at each iteration step 
from the minimum of the free energy. Since we neglect 
the Coulomb interaction in this part of the treatment, we 
can also eliminate the internal energy U. The free energy 
normalized as 



fs-U + /kin + /m 



i7T,2(l-t2)(i-i4)' 
in the dimensionless representation reads 



/ 



l-<2 

-(V X Q 



(l-i2)(l_t4) 



(89) 



When we substitute 
respect to c, d{f)/dc — 



4K^a; 

i^) into (p^, the variation with 
0, yields the condition for c as 

, (V^)'\ 



l-t2 



4^2^; 



1 - ^2 ^l-(l+c)L^(l-t4) / ■ 



(90) 



Condition (^) is not convenient for the numerical 
treatment. When the starting value of the wave func- 
tion is reasonable, or after a few iteration steps (a), (b), 
(a), (b), the correction c will be small, c ^ 1, so that 
the linear approximation of d90) is sufficient. 



2 (s - UJ ^ ujQ'^ ^ g) 



i2(l + f2) ^^2 (1 _ ^(1 _ i4))-iy 



(91) 



As a third iteration step (c) we thus may use Eq. ( p7| ) 
with c given by (|9l|). The iteration procedure we use 
starts from a preliminary adjustment of the wave function 
by a few steps, (a), (c), (a), (c), . . ., putting 5k = 0. 
After that the full iteration scheme (a), (c), (b), (a), (c), 
(b), ... is applied yielding all Fourier coefficients ok, 5k. 



E. Magnetic properties 



D. Accelerated iteration scheme 



As shown in Rcf. 34, the convergence is accelerated 



if the amplitude of the wave function is optimized after 



In Figs. 1-4 we present some numerical results to illus- 
trate the properties of the Bardeen equations. As one can 
see from the dimensionless equations (fz^ ) and (75), the 
behavior of the system is determined by a single mate- 
rial parameter, the GL parameter k. We assume Niobium 
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doped with non-magnetic impurities of a density giving 
the GL parameter kq = 1.5. 




-0.5 -0.5 



FIG. 1. The condensate density (plotted as uj{x,y)) in the 
triangular lattice for temperature t = 0.5, magnetic induction 
B /Bc2 =0.5, and GL parameter kq = 1-5. In the vortex cen- 
ters the condensate density ns{x,y)/n = (1 — t'^)Lo{x,y) goes 
to zero. Between the vortices uj{x,y) approaches its equilib- 
rium value 1 (which would be constant in the absence of a 
magnetic field) yielding nf^/n = 1 — t** = 0.94. 

Figure [l] shows a fishnet plot of the condensate density 
ns{x,y)/n = (1 — t^)uj{x,y) for temperature t — T/Tc = 
0.5 and the mean magnetic field B = 0.5 Bc2- The dips 
in Us correspond to the nodes of the wave function ip 
located at the vortex centers. The condensate density 
reaches its maximum between the vortices. 




-0.5 -0.5 



FIG. 2. The magnetic field in units of the upper critical 
field Bc2 for t = 0.5, 3/3^2 = 0.5, and ko = 1.5 as in Fig. 0. 
B(x,y) reaches its maximum Bmax at the vortex centers. 



Note that the condensate density is smaller than its 
non- magnetic value, = n(l — T'*/T^), also on the bor- 
ders of the elementary cells where the current is zero. 
This shows that non-local effects given by gradient cor- 
rections, e.g. the second term of (|63|), are important not 
only at the vortex core but also between the vortices. 

A complementary picture offers the plot of the mag- 
netic field B presented in Fig. The magnetic field 
reaches its maximum value, i?max, at the vortex cen- 
ters. This maximum field is very close to but slightly 
higher than the applied field Ba because the supercon- 
ductor tries to expel the magnetic field and compresses it 
into vortices. The magnetic pressure on the condensate 
is one of the forces balanced by the electric field. 




X, y 

FIG. 3. Profiles of the condensate density n{x, y) at various 
temperatures for B/Bc2 = 0.5 and = 1.5. The solid lines 
show cuts along the x-axis, and the dashed lines along the 
J/- axis. 



The temperature dependence of the condensate den- 
sity ns{x,y)/n = (1 — t'^)uj{x,y) is shown in Fig. |[ As 
one expects, the density of the condensate decreases as 
the temperature approaches the critical value, t ^ \. 
The dominant part of this decrease can be attributed to 
the reduced fraction of the condensate expressed by the 
factor (1 — i^). We have numerically checked that near 
the critical temperature the condensate density results 
identical to the solution of the standard GL theory. 

In Fig. ^ the density of the condensate is normalized to 
its value in absence of the magnetic field, ns(x,y)/n^ = 
u;{x,y). The suppression of uj{x^y) in the region be- 
tween the vortices is completely due to the magnetic field. 
One can see that at lower temperatures the condensate is 
less suppressed than it would result from the GL theory, 
where the latter one is equal to the curve at t = 0.99. 
This follows from the fact that the condensation energy 
Ts increases with the condensate density slower than the 
quadratic function of the GL theory. 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

X, y 

FIG. 4. Profiles of tlie reduced condensate density u){x,y) 
at various temperatures t = T/Tc for B / Bc2 = 0.5 and 
Ko = 1.5 as in Fig. ^ 



F. Virial theorem 

In the above treatment, the magnetic field was speci- 
fied by the mean value B of the magnetic induction in the 
sample. Macroscopic magnetic properties of the system, 
however, are given by the magnetization Af . Let us link 
these two quantities. 

For simplicity we assume that the sample is an infi- 
nite cylinder in the direction of the applied field. In this 
longitudinal geometry one has B = Ba + M. Since the 
non-local terms of the free energy are terminated at the 
second-order derivatives (the term called kinetic energy), 
we can conveniently use the mrial theorem derived by Do- 
ria, Gubernatis, and RaineiE3 and generalized by Klein 
and P6ttingerE3, to evaluate the applied magnetic field 

Ba- 

The idea of the virial theorem is as follows. Let us 
introduce a parameter i, which scales coordinates x and y. 
With this scaling one can generate a new wave function 
U!'{r) = u!{Lr). Since the mean magnetic field is given by 
the density of vortices, it scales as B'{r) — t~^i3(tr). We 
rescale all magnetic fields with except for the applied 
field Ba which is an external parameter. 

From B = V X A one can see that the vector potential 
scales as A'(r) = i~-'^A(tr), i.e., in the same way as a 
gradient. Accordingly, the density of kinetic energy scales 
with The mean value of free energy corresponding 
to the new wave function reads 



(/') 



-t- r 



UJ t^^Jl - L0{1 - t^) 



(92) 



The condensation energy (the first term) is independent 
of the scaling. The kinetic energy (the second term) 
scales with The magnetic energy (the third term) 
has three contributions: (6^) = ((V x Q)^) which scales 
with —2bba = — 2((V x Q)6a) which scales with t"^; 
and ba which is independent of the scaling. 

Since the scaling deforms the wave function and the 
internal magnetic field from their equilibrium values, the 
free energy (/') is greater than the free energy (/). For 
L = 1 the free energy (/') reaches its minimum being 
equal to (/). This minimum is given by a variation with 
respect to t, 



l-if) 



= 0. 



Condition 



2hbn 



in the explicit form, 



ujQ 



2((V X Q)' 



(93) 



(94) 



is called the virial theorem. Since 6, w and Q are known, 
the virial theorem ( p4[ ) provides us with the value of 
the applied magnetic field ba without having to take the 
derivative of the computed free energy. 

A convenient form of the virial theorem valid only 
for the Bardeen equations makes use of the Schrodinger 
equation (u4) from which follows 



s) 



(95) 



with s from (| 



The applied magnetic field then reads 
(262 + ^ _ s) 



26 



(96) 
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FIG. 5. The magnetization —M = Ba~Basa, function 
of the applied magnetic field Ba in units of the upper critical 
field Bc2 at temperatures t = 0.999, 0.85, 0.7, 0.5, 0.3 for 

Ko = 1.5. 
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The magnetization —M — Ba — B as a function of the 
applied magnetic field Ba is shown in Fig. || for different 
temperatures. At temperatures close to Tc the magneti- 
zation follows the line well known from the GL theory. 
Below the lower critical field B^i the system is in the 
Meissner state and —M — Ba- Above Bd the magneti- 
zation decreases and linearly vanishes at Bc2 where the 
system undergoes a transition into the normal state. 

At very low temperatures, the magnetization is de- 
formed into an S-shape. The slope of the decrease, 
dM/dBa, close to Bc2 increases with decreasing temper- 
ature and at a certain temperature Ta becomes infinite. 
Below Ta, the magnetic behavior of the system achieves 
an anomal feature. As the magnetic field is lowered from 
some high value, the system undergoes a first order tran- 
sition from zero to a finite magnetization at a field which 
is above Bc2- Since the free energy of the system with 
finite magnetization is lower than the free energy of the 
normal state, the system jumps to a finite magnetiza- 
tion as soon as the applied magnetic field allows for such 
solution. 

Such anomalous magnetic-transition has been observed 
by Ehrat and RindereiOH for Lead doped with Nio- 
bium. In spite of this experimental result we believe that 
the first order transition seen in Fig. ^ is an artifact of 
the Bardeep-.cypproximation. Indeed, detailed theoretical 
discussiong^acll of this anomalous behavior point to the 
important role of scattering on impurities. This mecha- 
nism is absent in the Bardeen approximation. 

The temperature Ta can be determined from the 
Bardeen equations. Close to the critical field Bc2 the 
density of condensate is small and one can expand the 
effective potential (|5^) into the GL form ( |5l| ) with coef- 
ficients 



2n^ 



Ti 



(97) 



For these asymptoticj^alues one can introduce an asymp- 
totic GL parameter,Ej 



(98) 



Let us return to features related to the electrostatic 
forces. As mentioned above, in the vortex core the mag- 
netic field B(x,y) is compressed and thus exceeds the 
value of the applied field Ba- In Fig. ^ we compare the 
applied field with the field Bmax in the center of the vor- 
tex. For all temperatures the compression is stronger at 
lower magnetic fields. 




As one can see from (|97|), this asymptotic GL parameter 
decreases with the temperature, 



FIG. 6. The applied field Ba (solid lines with dots) and 
the field Bmax in the vortex center (solid lines with crosses) 
plotted versus the induction B for ko = 1.5 at temperatures 
t = 0.999, 0.85, 0.7, 0.5, 0.3 as in Fig. |. For clarity, each line 
pair of the next temperature is shifted up by 0.2. The dashed 
lines indicate the large-K- limit Ba = B. 



VI. ELECTROSTATIC POTENTIAL AND 
CHARGE IN THE ABRIKOSOV LATTICE 

The electrostatic potential, (p, together with the 
Sorokin thermodynamic potential, Ws, control the mo- 
tion of Cooper pairs. Indeed, the total effective potential 
acting on the Cooper pairs is x = e*V + ^Ws- The sep- 
aration of the effective potential x into its electrostatic 
and thermodynamic components sheds a light on the role 
of the electrostatic potential in the Schrodinger equation 
or {M). 



Kot- 



(99) 



The transition temperature Ta appears when the effective 
GL parameter equals to 1/^/2, i.e.. 



T, 



T, 



V2, 



Ho 



(100) 



For Ko — 1.5 one finds Ta = 0.47 Tc. We expect that one 
should be cautious about results of the Bardeen equations 
below Ta- 



A. Electrostatic potential 

The electrostatic potential in the vortex lattice is given 
by the screened Poisson equation (p9). For simplicity 
we neglect the screening, putting A^\7^eip ~ 0. This 
approximation is justified below. 

To be compatible with the above notation, we define a 
dimensionless electrostatic potential, 



en 



(101) 
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With the screening neglected one finds from (p9| ) 

(j) = s - Lu + CiLO + C2\/l - (l - t^)uj, (102) 



with temperature dependent factors 

1 91n£con 



and 



C2 = 



1 — <2 9 Inn 
1 9 In 7 



d\ni 



(103) 



(104) 



The term s — a; in (102) corresponds to xlV'l /^^ in (|59|). 
The Ci^2 terms correspond to the second and third terms 
of (p9|), respectively. 




-0.5 -0.5 



FIG. 7. The electrostatic potential 4>{x,y), Eq. (101). The 
temperature t — 0.5, the magnetic field B / Bc2 ~ 0.5 and the 
GL parameter kq = 1-5 are identical to the values used in 
Figs. and|^. The thermodynamic coefiicients Ci and C2 are 
specified in Tab. 1. 

Figure ^ shows a fishnet plot of the electrostatic po- 
tential. The potential reaches its minimum at the vortex 
centers, i.e., it attracts electrons to vortices. 

The total electrostatic potential is composed of three 
components: the Bernoulli potential (f>B = s — w, the 
contribution due to the condensation energy = Ciw 
and the reduced normal state thermoelectric potential 
02 = C2-\/l — (1 — t^)w. Individual components are 
compared in Fig. ^. 

The Bernoulli potential 0^ shown in Fig. || is negative. 
Due to the quasiparticle screening, the Bernoulli poten- 
tial reaches zero at the center of the vortex. With respect 
to the center of the vortex, the forces corresponding to 
the Bernoulli potential are repulsive inside the core while 
they are attractive outside. 

The potential 01 caused by the density dependence of 
the condensation energy is positive. Being proportional 



to the density of condensate, it has a minimum at the 
vortex center where it reaches zero. For Niobium this 
contribution is dominant since the coefficient Ci = 1.9 is 
rather large compared to the coefficients of other contri- 
butions. 




0.4 0.6 
x/a, y/a 

FIG. 8. The components of the electrostatic p oten tial in 
the vortex lattice of spacing a according to Eq. (lOS). The 
individual potentials are: the total potential </!> (solid lines), 
the Bernoulli potential (j>B = s — u; (dashed lines), the conden- 
sation potential <^i — C\lj (dashed-dotted lines), and the nor- 
mal thermodynamic potential (/!>2 ~ C2\/l — (1 — t'*)aj (dot- 
ted lines). The splitting of lines at larger distances charac- 
terizes the x-direction (lower curves) or y-direction (upper 
curves). Parameters as in Fig. 7. 

The normal state thermodynamic potential 02 is also 
positive giving the only non-zero contribution at the vor- 
tex center. One can see that 02 reduces the total po- 
tential since it has the maximum at the vortex core and 
falls outside. Its coefficient C2 = 0.42 is about four times 
smaller than Ci, therefore this term cannot cancel the 
potential 0i. 

We want to stress that even at temperature t = 0.5 
when 96% of electrons are in the condensate, the thermo- 
dynamic correction to the electrostatic potential cannot 
be neglected. This result contradicts the temperature 
dependence, of the thermodynamic correction derived by 
Rickayzencj, see Eq. (|3|) . Within the hydrodynamic pic- 
ture one can show that the limit of weak currents adopted 
by Rickayzen is responsible for this disagreement. In this 
limit, the effect of the diamagnetic current on the con- 
densate density rig vanishes as T — ^ 0, see Eq. (68). The 
temperature dependence of the thermodynamic correc- 
tion merely reflects the temperature dependence of n^. 
The limit of weak current does not apply to the vortex 
core. In the center of the vortex core, the condensate 
density has to go to zero keeping the magnitude of the 
thermodynamic correction appreciable at any tempera- 
ture. 
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B. Effective potential 



C3 = ^^(l-t2)(l-i4). 



(107) 



To enlighten the role of the electrostatic potential in 
the balance of forces in superconductors, we compare the 
effective potential the electrostatic potential acting on 
the Cooper pair e*(^, and the thermodynamic potential 
of Sorokin 210,^ in Fig. ^ All these contributions are in 



dimensionless units corresponding to (101). 

One can see that the electrostatic potential is not a 
small correction to the effective potential of a thermo- 
dynamic origin. The amplitude of the electrostatic po- 
tential is about an order of magnitude larger than the 
amplitude of the effective potential x- Accordingly, the 
effective potential x = e.*^ + 2u;s results from a strong 
compensation of the thermodynamic potential 2ws and 
the electrostatic potential e*(p. 




0.4 0.6 

x/a, y/a 

FIG. 9. The effective potential x ~ e'ip + 2vjs (full 
lines), the electrostatic potential acting on the Cooper pair 
e*95 (dashed lines) and the thermodynamic potential 2ws 
(dash-dotted lines) for kq = 1-5. Parameters and presentation 
as in Fig. H 



C. Charge 



The distribution of the charge in the vortex lattice is 



given by the Poisson equation, p - 
a dimensionless charge. 



P_ 

en ' 



-eV iy9. We introduce 



(105) 



which measures the relative deviation of the charge den- 
sity from the crystal value. In dimensionless representa- 
tion the Poisson equation reads 



X2 



•^Lon 



(106) 
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FIG. 10. The function —V^tf) proportional to the charge 
density p{x,y). The amplitude of the dimensionless charge 



density is p = 
as in Figs. 7, 



9.510' 
8, 9. 



■{\—t )(1— t ) V 4>- Same parameters 



Figure |To| shows a fishnet plot of the charge distri- 
bution. In the vortex core the charge is depleted, the 
missing charge is distributed between vortices. 

A striking feature is the very rapid change of the charge 
sign at the distance about 0.4 from the vortex center. 
While the charge in the core is rather flat, its spatial 
variation between vortices is quite strong. This picture 
of the charge distribution is just opposite-|to the one as- 
sumed by Kumagai, Nozaki and MatsudaE2l who expected 
a flat charge distribution between vortices. Comparing 
these two pictures, however, one has to keep in mind that 
Kumagai et al discuss YBCO with k, ~ 100 while Fig. |l0| 
presents the case of kq = 1.5 in Niobium. 

The particular shape of the charge seen in Fig. ^ re- 
sults from the interplay between the Bernoulli potential 
(j)B and the potential 0i due to the condensation en- 
ergy. In Fig. |11|, we show the charge density decomposed 
into contributions corresponding to individual potentials, 

Pi CX V^0j. 

The charge distribution pB corresponding to the 
Bernoulli potential (pB has its maximum at the center of 
the vortex. In Ref. |2^, where only the non-local (quan- 
tum) Bernoulli potential has been assumed, there is a 
minimum of the charge density in the center of vortex. 
The local maximum seen in Fig. 



11 



with 



follows from the 
quasiparticle screening not assumed in Ref. ^ 

The amplitudes of the contributions pi and p2 depend 
on the constants Ci and C2, which strongly depend on 
the material in question. For Niobium one has Ci = 1.9 
and C = 0.42; therefore pi dominates. In Appendix A 
one can see that both constants Ci and C2 are propor- 



15 



tional to the slope of the density of states at the Fermi 
level. In general one can say that the amplitude of p2 is 
smaller than the amplitude of pi and the two contribu- 
tions have opposite signs of the periodic parts. We note 
that due to the large value of Ci and small C2 , the total 
charge has a minimum in the center of the vortex. 




-lOoL:^ , , , , liji 

0.2 0.4 0.6 0.8 1 
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FIG. 11. The components of the charge density in the 
notation of Fig. ^ for same parameters. The individual 
charge densities are: the total p (solid lines), the Bernoulli pB 
(dashed lines), the condensation part pi (dash dotted lines), 
and the normal thermodynamic part p2 (dotted lines). 



D. Screening and the quasi- neutral approximation 

For Niobium, the Thomas-Fermi screening length is 
very small, 

^=2.510-6(1-^4). (108) 

One can thus neglect the screening, ^TF^hon^^'l^ ^ ^■ 
Indeed, the Laplace operator in the Fourier representa- 
tion is ^ if^ = (47r/\/3)K6(«2 + ij Since k is 
of the order of unity, 6 < 1 and the number of needed 
Fourier components is also limited, i,j < 100, the screen- 
ing is negligible for all Fourier components considered. 
For Niobium, the factor 

^^^ = 3.810-5 (109) 

which determines C3 in (|l06| ), is also very small. Sim- 
ilarly small value can be expected for any conventional 
superconductor. It leads to relative charges of the order 
of 10^^°. The quasi-neutral approximation, j{n) « ^(no) 
etc., is thus well justified when one solves for the wave 
function and the vector potential. 



VII. CONCLUSIONS 

In this paper we have discussed the electrostatic po- 
tential in the Abrikosov lattice of vortices. To this end 
we have derived a set of three Ginzburg-Landau equa- 
tions which include the Maxwell equation for the vector 
potential, the Schrodinger equation for the wave func- 
tion, and the Poisson equation for the electrostatic po- 
tential. These equations determine the minimum of the 
free energy made of four components: the condensation 
energy of Gorter and Casmir; the quantum kinetic en- 
ergy of Ginzburg and Landau; the magnetic free energy 
of Helmholtz; and the Coulomb energy. 

The marriage of the Gorter-Casimir two-fluid model 
with the Ginzburg-Landau theory has been suggested 
earlier by Bardeen who has also discussed properties of 
this theory at different temperature limits. We have em- 
ployed his approach as it offers a very simple extension 
of the GL theory towards low temperatures. As our re- 
sults document, this extended theory can be treated with 
standard numerical tools of the GL theory. 

With the electrostatic interaction included, the effec- 
tive potential acting on the superconducting condensate 
is naturally a sum of the electrostatic potential and the 
thermodynamic potential. One can say that the electro- 
static potential (over-)screens the thermodynamic poten- 
tial leaving a relatively small effective potential. 

In spite of the very important role of the electrostatic 
potential among forces acting on the condensate, the elec- 
trostatic potential can be eliminated from the Ginzburg- 
Landau theory so that one has to solve a set of two, not 
three, equations. This simplification is possible by two 
reasons. First, the charge modulation which corresponds 
to this potential, is so small on the scale of the charge 
density in metals that one can neglect its effect on local 
values of material parameters. 

The second reason is more fundamental. As noticed 
by van Vijfeijken and Staas, there is a force between the 
condensate and the normal electrons. This force keeps 
the normal electrons at rest, i.e., it balances the electric 
field having an equal amplitude and the opposite orienta- 
tion. The force of van Vijfeijken and Staas is an exclusive 
function of the condensate density. Accordingly, one can 
express the electric force or the electrostatic potential as 
a function of the condensate density. In this way the 
electrostatic potential can be unified with the thermody- 
namic potential into an effective potential of GL-type. 

In the numerical treatment we have used the param- 
eters of Niobium. Our choice of this conventional ma- 
terial was determined by known empirical rules needed 
to predict amplitudes of the individual contributions to 
the electrostatic potential. We expect that other d-band 
superconductors behave similarly. 

Finally, we would like to stress that the presented the- 
ory is simplified in many directions. First, it is restricted 
to isotropic materials. We have omitted all features of 
the band structure except for the density of states and 
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its slope on the Fermi level. Second, the two-fluid model 
of Gorter and Casimir describes only gross features of the 
thermodynamics of superconductors. Third, the gradient 
approximation of the Ginzburg-Landau theory is justified 
only close to the critical temperature, at low tempera- 
tures one has to take the kinetic energy of Ginzburg and 
Landau as an ad hoc approximation. In the future, we 
plan to address layered structures and use a more general 
form of the Gorter-Casimir model. 



Now we can complete the estimate of d^/dn. From 



(A.3 A5) follows 



dV 



= (1 + 2A) 



therefore relation (A2) can be expressed as 

2A91nX'o 



dn 



1 + X dEp 



(A6) 



(A7) 



APPENDIX A: ESTIMATE OF MATERIAL 
PARAMETERS FOR NIOBIUM 



2. Coefficient decon/dr, 



In this appendix we estimate material parameters, 
d^/dn and dscon/dn, which determine the electrostatic 
potential in the superconductor, see ( p9| ) . To be specific 
we assume Niobium. 



1. Coefficient dj/dn 

The linear coefficient of the specific heat 7 is linked to 
the density of states V per spin and unitary volume. 



7 = g^^^B 



V. 



(Al) 



It is advantageous to express the density derivative of 7 
in terms of the energy derivative of the density of states. 
Using dEp/dn = 1/2V we find 



97^1^2fc.ainP 



dn 



OEf 



(A2) 



The density of states T) includes the mass mnormaliza- 
tion due to the electron-phonon interaction,c3 



(A3) 



where Vq is a bare density of states and A is the coupling 
parameter. The value and the energy derivative of Vq is 
provided by ab initio studies of Niobium.E^ 

The value of the coupling parameter A is found com- 
paring V from the experimental 7 with the theoretical 
Vq. The energy derivative of A, however, is not provided 
in the literature. To estimate the derivative of A we write 
it as a product, 



A = VoV, 



(A4) 



where V is the BCS interaction. 

According to trends found from the effects of impuri- 
ties on the critical temperature and the specific heat, the 
major changes of A follow from the density of states, while 
the BCS interaction V remains nearly constant .E3 As a 
first approximation we thus assume 



dV 
dn 



= 



dV 
dEp 



= 0. 



(A5) 



The derivative of the condensation energy (^0|) includes 
the derivative of the critical temperature. For Niobium 
and similar materials the critical temperature is given by 
the McMillan formula,El 



T, 



1.45 



exp 



1.04(1 -f A) 
A-M*(l + 0.62A) 



(A8) 



where 9d is the Debye temperature and fi* is the 
Coulomb pseudopotential. From ( ^ ) and (A8) we ex- 
press the condensation energy as 



12.6 



fc|(l + A)X>o6'^ exp 



-2 



1.04(1 + A) 
A-/i*(l + 0.62A) 

(A9) 



Experience from dilutej-aJloys shows that the product 
VqO^ is nearly constant.Q We thus use as the second 
approximation, 



0. 



(AlO) 



In this approximation the derivative of the condensation 
energy is given by the derivative of the factor 1 + A and 
by the derivative of the argument of the exponential. 



dec 



dn 



d_ 

' dn 



2.08(1 + A) 
A-/i*(l + 0.62A) 



ln(l 



(All) 



Again, the experience with dilute alloys shossfs that the 
Coulomb pseudopotential is nearly constant jlfl therefore 
we take as the third approximation. 



d^jL* 
dn 



(A12) 



With approximation (A12) the density derivative of the 



condensation energy becomes proportional to the deriva- 
tive of the coupling parameter. 



de. 



dn 




(A13) 
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The den sity derivative of the cou phng constant follows 
from (A4) and approximation (A5) as 



dX 
dn 



V dlnVo 
2(1 + A) dEp 



(A14) 



The derivative of the condensation energy is thus propor- 
tional to the BCS interaction, 



dec 



EconV dlnVo / 1.04(l + 0.38/r)(l + A) 1 



dn 



(1 + A)2 dEp \^ (A-/x*(l + 0.62A))' 2j' 

(A15) 



The material parameters for Niobium which we have 
used are listed in Table 1. For convenience, we have 
included values which can be evaluated from the above 
formulas, e.g., the critical temperature is given by (|A8|). 
The logarithmic derivative of the density of states with 
respect to the energy is extracted from the figure in 
Ref. The hole density n_has been evaluated from 
the London penetration depthc3 



A 



(A16) 



At zero temperature all holes are in the condensate, n = 
Us- The listed density of holes follows from Aloh = Aq = 
3.9 10~^ m and the mass mo = 1.2 mg. This effective 
mass is an estimate of values 1.12, l£j 1.28 and 1.22 for 
different orbits of the pure Niobium. CJ 

We assume that the properties of the material are mod- 
ified by oxygen impurities of a concentration ranging 
from to 0.03. We neglect the effect of impurities on 
the thermodynamic parameters taking into account only 
their dominant effect on the London penetration depth 
and the GL coherence length. In the dirty limit, the GL 
coherence length, defined in our model as 



e = 



(A17) 



scales with the square- root of the mean free path ^ oc 
a/I, while the effective London penetration depth scales 
with its inverse, Aloh oc 1/-\/7E3 Accordingly, the GL 
parameter n = Aloh/C is proportional to the inverse mean 
free path, k oc l/l. One can see that the proper scaling 
of both characteristic lengths is achieved by the scaling 
of the effective mass, 



m — niQ- 



^pure 



(A18) 



where Kpure is the GL parameter of the pure Niobium 
while Ko is the actual value for a given concentration of 
impurities provided in Ref. [42. 



critical temperature''^ 

Debye temperature"*^ 9o 

coupling parameter^^ A 

Coulomb pseudopot.'*^ fi* 

coef. of spec, heat''^ 7 

mass in pure Nb^^ mo 

hole density'*^ n 
log. der.^'^ 

GL parameter^^ Kpuio 

density of states (Al) V 5.7 lO'^'^ J^^m^^ 

bare density . . . (A3) Vq 3.0 10"*^ J^^m-^ 

BCS interaction (A4) V 2.9 lO^^** Jm^ 

cond. energy (30) Scon 1-6 iC Jm~^ 

cond. en. per pair 9.17 10^^ eV 

coefficient (A7) ifjT'c 3.85 10"^ eV 



9.5 K 
275 K 
0.89 
0.15 

719 Jm-3K-2 

1.2 TTie 

2.2 10^8 m-3 



1.1 1019 j-i 



0.78 



coefficient (A15) 
coefficient of Ci 



%^ 8.73 10-*^ eV 

on 

1.9 

In n 



coefficient of C2 0.42 



Table 1. Material parameters of pure Niobium. 
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